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Abstract. We examine some numerical iterative methods for computing the eigenvalues and eigenvectors 
of real matrices. The five methods examined here range from the simple power iteration method to the more 
complicated QR iteration method. The derivations, procedure, and advantages of each method are briefly 
discussed. 



1. Introduction 

Eigenvalues and eigenvectors play an important part in the applications of linear algebra. The naive 
method of finding the eigenvalues of a matrix involves finding the roots of the characteristic polynomial of 
the matrix. In industrial sized matrices, however, this method is not feasible, and the eigenvalues must 
be obtained by other means. Fortunately, there exist several other techniques for finding eigenvalues and 
eigenvectors of a matrix, some of which fall under the realm of iterative methods. These methods work by 
repeatedly refining approximations to the eigenvectors or eigenvalues, and can be terminated whenever the 
approximations reach a suitable degree of accuracy. Iterative methods form the basis of much of modern 
day eigenvalue computation. 

In this paper, we outline five such iterative methods, and summarize their derivations, procedures, and 
advantages. The methods to be examined are the power iteration method, the shifted inverse iteration 
method, the Rayleigh quotient method, the simultaneous iteration method, and the QR method. This paper 
is meant to be a survey over existing algorithms for the eigenvalue computation problem. 

Section 2 of this paper provides a brief review of some of the linear algebra background required to 
understand the concepts that are discussed. In section 3, the iterative methods are each presented, in order 
of complexity, and are studied in brief detail. Finally, in section 4, we provide some concluding remarks and 
mention some of the additional algorithm refinements that are used in practice. 

For the purposes of this paper, we restrict our attention to real-valued, square matrices with a full set of 
real eigenvalues. 

2. Linear Algebra Review 

We begin by reviewing some basic definitions from linear algebra. It is assumed that the reader is 
comfortable with the notions of matrix and vector multiplication. 

Definition 2.1. Let A G M™ xn . A nonzero vector x G R™ is called an eigenvector of A with corresponding 
eigenvalue A G C if Ax — Ax. 

Note that eigenvectors of a matrix are precisely the vectors in R™ whose direction is preserved when 
multiplied with the matrix. Although eigenvalues may be not be real in general, we will focus on matrices 
whose eigenvalues are all real numbers. This is true in particular if the matrix is symmetric; some of the 
methods we detail below only work for symmetric matrices. 

It is often necessary to compute the eigenvalues of a matrix. The most immediate method for doing so 
involves finding the roots of characteristic polynomials. 

Definition 2.2. The characteristic polynomial of A, denoted Pa{%) for x G R, is the degree n polynomial 
defined by 

P A {x) := det(x/ - A). 
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It is straightforward to see that the roots of the characteristic polynomial of a matrix are exactly the 
eigenvalues of the matrix, since the matrix zl — A is singular precisely when z is an eigenvalue of A. It 
follows that computation of eigenvalues can be reduced to finding the roots of polynomials. Unfortunately, 
solving polynomials is generally a difficult problem, as there is no closed formula for solving polynomial 
equations of degree 5 or higher. The only way to proceed is to employ numerical techniques to solve these 
equations. 

We have just seen that eigenvalues may be found by solving polynomial equations. The converse is also 
true. Given any monic polynomial 

f(z) = z n + a n _iz™ _1 + . . . + a x z + a , 

we can construct the companion matrix 

" 
1 
1 

1 -a„-2 
1 -a„_i 

It can be seen that the characteristic polynomial for the companion matrix is exactly the polynomial 
f(z). Thus the problem of computing the roots of a polynomial equation reduces to finding the eigenvalues 
of a corresponding matrix. Since polynomials in general cannot be solved exactly, it follows that there is no 
method that will produce exact eigenvalues for a general matrix. 

However, there do exist methods for computing eigenvalues and eigenvectors that do not rely upon solving 
the characteristic polynomial. In this paper, we look at some iterative techniques used for tackling this 
problem. These are methods that, when given some initial approximations, produce sequences of scalars or 
vectors that converge towards the desired eigenvalues or eigenvectors. On the other hand, we can make the 
notion of convergence of matrices precise as follows. 

Definition 2.3. Let A^\ A^ 2 \ A^\ . . . be a sequence of matrices in R" iX ™. We say that the sequence of 
matrices converges to a matrix A G u mxn if the sequence A+j of real numbers converges to Aij for every 
pair 1 < i < m, 1 < j < n, as k approaches infinity. That is, a sequence of matrices converges if the 
sequences given by each entry of the matrix all converge. 

Later in this paper, it will be necessary to use what is known as the QR decomposition of a matrix. 

Definition 2.4. The QR decomposition of a matrix A is the representation of A as a product 

A = QR, 

where Q is an orthogonal matrix and R is an upper triangular matrix with positive diagonal entries. 

Recall that an orthogonal matrix U satisfies U T U = /. Importantly, the columns of Q are orthogonal 
vectors, and span the same space as the columns of A. It is a fact that any matrix A has a QR decomposition 
A = QR, which is unique when A has full rank. 

Geometrically, the QR factorization means that if the columns of A form the basis of a vector space, then 
there is an orthonormal basis for that vector space. This orthonormal basis would would form the columns 
of Q, and the conversion matrix for this change of basis is the upper triangular matrix R. The methods for 
obtaining a QR decomposition of a matrix has been well studied and is a computationally feasible taskQ. 

At this point, we turn our attention to the iterative methods themselves. 



i The simplest method for computing a QR factorization of a matrix A is to apply the Gram-Schmidt algorithm on the 
columns of A. A student of linear algebra may be horrified at the prospect of carrying out this tedious procedure once, let 
alone once per iteration of an iterative method, but recall that when working with matrices of thousands of rows or more, 
all computations are done electronically. Furthermore, the methods used to calculate a QR decomposition are usually more 
complicated than Gram-Schmidt; for instance, a common technique is to apply a series of Householder transformations [2j. 
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3. Description of the Iterative Methods 



The iterative methods in this section work by repeatedly refining estimates of eigenvalues of a matrix, 
using a function called the Rayleigh quotient. 



Definition 3.1. Let A S R nx ™. Then the Rayleigh quotient of a nonzero vector x £ R n is 

x 1 



" T Ax 



r(x) 

X' X 

Note that if x is an eigenvector for A with corresponding eigenvalue A, then the Rayleigh quotient for x is 

. . x A.X Xx x 

X X X X 

which is exactly the corresponding eigenvalue. 

In fact, given any nonzero x G R™, the Rayleigh quotient r(x) is the value that minimizes the function 
f(a) = \\ax — Ax\\ 2 over all real numbers a, which measures the error incurred if x is assumed to be an 
eigenvector of A. Thus, if we treat x as an estimate for an eigenvector of A, then r(x) can be seen as the 
best estimate for the corresponding eigenvalue of A, since it minimizes this error value. 

We are now ready to consider the first technique for iterative eigenvalue computation. 

3.1. Power Iteration. Let A £ R nx ™. Recall that if q is an eigenvector for A with eigenvalue A, then 
Aq = Xq, and in general, A k q = \ k q for all k € N. This observation is the foundation of the power iteration 
method. 

Suppose that the set {qi} of unit eigenvectors of A forms a basis of R™, and has corresponding real 
eigenvalues {A^} such that |Ai| > |A 2 | > . . . > |A„|. Let be an approximation to an eigenvector of A, 
with ||y(°)|| = 1. (We use the term "approximation" in this situation quite loosely, allowing it to refer to 
any quantity that is believed to be "reasonably close" to the correct valu^l.)We can write as a linear 
combination of the eigenvectors of A] for some ci, . . . , c„ G R we have that 



cxqi + . . . + c n q ni 



and we will assume for now that ci =/= 0. 
Now 



Aw (0) = ciAigi + c 2 A 2 <?2 + ■ • ■ + c n X n q n 



and so 

A k v (o) 



Since the eigenvalues are assumed to be real, distinct, and ordered by decreasing magnitude, it follows 
that for alH = 2, . . . , n, 

K ^ k 




lim -i =0. 

fc— >oo \M J 

So, as k increases, A k v^ approaches c\X k qi, and thus for large values of k, 

A k v^ 

The method of power iteration can then be stated as follows: 

Pick a starting vector with = 1 

For k = 1,2, ... 
Let w = Av^ 1 " 1 



2 These approximations are always made when the correct value is unknown (indeed, they are made in an attempt to determine 
the correct value); we therefore do not require any particularly demanding conditions on how "close" the approximation must 
be to being correct. 



Let i/ fc) = w/ \\w\\ 

In each iteration, gets closer and closer to the eigenvector q±. The algorithm may be terminated 
at any point with a reasonable approximation to the eigenvector; the eigenvalue estimate can be found by 
applying the Rayleigh quotient to the resulting v^ k \ 

The power iteration method is simple and elegant, but suffers some major drawbacks. The method only 
returns a single eigenvector estimate, and it is always the one corresponding to the eigenvalue of largest 
magnitude. In addition, convergence is only guaranteed if the eigenvalues are distinct — in particular, the 
two eigenvalues of largest absolute value must have distinct magnitudes. The rate of convergence primarily 
depends upon the ratio of these magnitudes, so if the two largest eigenvalues have similar sizes, then the 
convergence will be slow. 

In spite of its drawbacks, the power method is still used in some applications, since it works well on large, 
sparse matrices when only a single eigenvector is needed. However, there are other methods that overcome 
the difficulties of the power iteration method. 

3.2. Inverse Iteration. The inverse iteration method is a natural generalization of the power iteration 
method. 

If A is an invertible matrix with real, nonzero eigenvalues {Ai, . . . , A„}, then the eigenvalues of A -1 are 
{1/Ai, . . . , 1/A„}. Thus if |Ai| > |A 2 | > . . . > |A„|, then |l/Ai| < |1/A 2 | < . . . < |1/A„|, and so by applying 
the power method iteration on A -1 , wc can obtain the eigenvector q n and eigenvalue A„. 

This gives a way to find the eigenvalue of smallest magnitude, assuming that A^ 1 is known. In general, 
though, the inverse matrix is not given, and calculating it is a computationally expensive operation. How- 
ever, computing x = A~ x b is equivalent to solving the system Ax = b for x given 6, and this operation can 
be efficiently performed. Fortunately, this is all that is required for the inverse iteration method, which we 
can now state as follows: 

Pick a starting vector with ||w ( - ' ) || =1 
For k = 1,2,... 

Solve Aw = for w 

Let = wj \w\ 

The advantage of inverse iteration is that it can be easily adapted to find any eigenvalue of the matrix 
A, instead of just the extreme ones. Observe that for any \i e R, the matrix B = A — [il has eigenvalues 
{Ai — fj,,...,X n — ^}. In particular, by choosing \i to be close to an eigenvalue Xj of A, we can ensure 
that Xj — yu is the eigenvalue of B of smallest magnitude. Then by applying inverse iteration on B, an 
approximation to qj and Xj can be obtained. 

This version of the algorithm, known as inverse iteration with shift, can be summarized as follows: 

Pick some fi close to the desired eigenvalue 
Pick a starting vector with =1 
For k = 1,2,... 

Solve (A - fj,I)w = v^ k ~^ for w 

Let v( k) = wj \w\ 

The inverse iteration with shift method allows the computation of any eigenvalue of the matrix. However, 
in order to compute a particular eigenvalue, you must have some initial approximation of it to start the 
iteration. In cases where an eigenvalue estimate is given, the inverse iteration with shift method works well. 

3.3. Rayleigh Quotient Iteration. The inverse iteration method can be improved if we drop the restriction 
that the shift value remains constant throughout the iterations. 

Each iteration of the shifted inverse iteration method returns an approximate eigenvector, given an esti- 
mate of an eigenvalue. The Rayleigh quotient, on the other hand, produces an approximate eigenvalue when 
given an estimated eigenvector. By combining these two operations, we get a new variation of the inverse 
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shift algorithm, where the shift value is recomputed during each iteration to become the Rayleigh quotient 
of the current eigenvector estimate. 

This method, called the Rayleigh quotient iteration method (or simply the RQI method), is as follows: 



Pick a starting vector v' ' with = 1 

Let AW =r(«<°>) := (f (0) ) T A(v^) 
For k = 1,2, ... 

Solve (A - X^-^^w = v^- 1 *) for w 

Let «W = w/ \\w\\ 

Let \W = r(v^) := (v^) T A(v^) 



In this method, we no longer need to have an initial eigenvalue estimate supplied; all that is required is 
an initial vector v^°\ The eigenvector produced depends on the initial vector chosen. Note that since each 
vector is a unit vector, we have (V^) — 1, simplifying the expression for r (V^). 

The main advantage of the RQI method is that it converges to an eigenvector very quickly, since the 
approximations to both the eigenvalue and the eigenvector are improved during each iteration. Thus far, 
we have not given any quantitative discussion on the speed at which an iteration sequence converges to the 
limit, since these technical details are not the purpose of this paper. However, it is worth mentioning here 
that the convergence rate of the RQI method is said to be cubic, which means that the number of correct 
digits in the approximation triples during each iteration [3J. In contrast, the other algorithms described in 
this paper all have the much slower linear rate of convergence^ 

One very significant disadvantage for the RQI method is that it does not always work in the general case. 
The method is only guaranteed to converge when the matrix A is both real and symmetric, and is known to 
fail in the cases where the matrix is not symmetric. 

3.4. Simultaneous Iteration. The methods discussed so far are only capable of computing a single eigen- 
value at a time. In order to compute different eigenvalues of a matrix, the methods must be reapplied 
several times, each time with different initial conditions. We now describe a method that is capable, in some 
situations, of producing all of the eigenvalues of a matrix at once. 

Once again, our starting point is the basic power iteration method. Let A be a real, symmetric, full rank 
matrix; in particular, A has real eigenvalues and a complete set of orthogonal eigenvectors. Recall that given 
the starting vector v^ ', we can write it as a linear combination of eigenvectors {qi} of A: 

v {0) = CiQi + • • • + c n q n . 

Note, however, that the power iteration method only looks at eigenvectors that are nontrivial components 
in this linear combination; that is, only the eigenvectors that are not orthogonal to have a chance of 
being found by the power iteration method. This suggests that by applying the power iteration method to 
several different starting vectors, each orthogonal to all of the others, there is a possibility of finding different 
eigenvalues. 

With this idea in mind, we may take the following approach. Begin with a basis of n linearly independent 
vectors {wj '* . . . } of R™ , arranged in the matrix 

y(0) = L(°) 

Let 



v (o) 



y{k) _ ^fcy(O) 



(fc) 



which effectively applies the power iteration method to all of the vectors {v^ . . . ul '} at once. We thus 
expect that as k — > oo, the columns of V^ k ' become scaled copies of qi, the unit eigenvector corresponding 
to the eigenvalue of largest magnitude. (Note that the columns are not necessarily unit vectors themselves, 
since we did not normalize the vectors in this version of the algorithm). 



^In further contrast, the well known Newton's method for rapidly finding roots of differentiable functions has quadratic 
convergence. Generally, the study of algorithm convergence rates is extremely technical, and the author feels no regret in 
omitting further details from this paper. 



So far, we have not found anything useful or new. We have obtained n eigenvectors, but they are possibly 
all in the same direction. The main development occurs when we decide to orthonormalize the columns 
of at each iteration. In the original power iteration method, the obtained eigenvector estimate was 
normalized in each iteration. In this multivector version, the analogue is to obtain an orthonormal set of 
eigenvector estimates during each iteration, forcing the eigenvector approximations to be orthogonal at all 
times. This is done by using the QR decomposition of V ^ . 

In each step of this iteration method, we obtain a new matrix, W, by multiplying A by the current eigen- 
vector approximation matrix, can then extract the orthonormal column vectors of Q from the 
QR decomposition of W, thus ensuring that the eigenvector approximations remain an orthogonal basis of 
unit vectors. This process is then repeated as desired. The algorithm, known as the simultaneous iteration 
method [3], can be written as follows: 

Pick a starting basis ,...,«„ } of M" . 



,(<>) 



„(0) 



Build the matrix V 

Obtain the factors 
For k = 1,2,... 
Let W = AQ( k -V 

Obtain the factors QWrW = W 

It is a fact that if the matrix A has n orthogonal eigenvectors q\ , . . . , q n with corresponding real eigenvalues 
| Ai| > . . . > | A n |, and if the leading principal submatrices of the product [qi\ ■■ ■ \q n ] T are nonsingular, 
then the columns of will converge towards a basis of eigenvectors of A. (Recall that the leading principal 
submatrices of the matrix B are the top left square submatrices of B.) So, at last, here is a method that, 
under some hypotheses, computes all of the eigenvectors of the matrix A at once. 

However, this method is generally not used in practice. It turns out that there is a more elegant form of 
this algorithm, which we examine now. 



3.5. The QR Method. The QR method for computing eigenvalues and eigenvectors pQ, like the simulta- 
neous iteration method, allows the computation of all eigenvalues and eigenvectors of a real, symmetric, full 
rank matrix at once. Based upon the matrix decomposition from which it earned its name, the simplest 
form of the QR iteration algorithm may be written as follows: 

Let A(°) = A 
For k = 1,2, ... 

Obtain the factors QWrW = A^ 1 ) 

Let AW =i?( fe )QW 

Simply put, in each iteration, we take the QR decomposition of the current matrix A^ k ~ x \ and multiply 
the factors Q and R in the reverse order to obtain the new matrix A^ . 

It is surprising that this non-intuitive process would converge to anything useful, let alone a full set of 
eigenvectors and eigenvalues of A. However, it turns out that this algorithm can, in some sense, be seen 
to be equivalent to the simultaneous iteration method. The simultaneous iteration method of the previous 
section can be written as follows: 

Let Q (0) = I 
For k = 1,2, ... 

Let W = AQ (k - 1) 

Obtain the factors 

Let AW = (Q( fe )) T AQW 
Let RW =RWR(k-i)... R W 



Here we have renamed the matrices of the previous section as Q , and have introduced the matrices 

and which do not affect the correctness of the algorithm. 
The QR method of this section can be rewritten in the following way: 

Let A® = A 
For k = 1,2,... 

Obtain the factors QWrW = A^" 1 ) 

Let AW =RWQW 

Let qM=QWqW-.-QW 

Let R {k) =i?( fc )i?( fe - 1 )..-i?( 1 ) 

Again, the introduction of the matrices Q^ k ' and R} k ^ do not affect the outcome of the algorithm. Note 
that by this algorithm, the following identities hold for all values of k: 

. A« = rWqW 
• (Q (fc) ) T Q (fe) =/ 

The similarity between these two algorithms becomes apparent by the following theorem. 

Theorem 3.2. The simultaneous iteration method and the QR method both generate the same sequences of 
matrices A^ k \ Q^ k \ and R} k \ which satisfy the following relations: 

(1) A k = Q^rW 

(2) A™ = (Q W ) T AQW 

Proof. We proceed by induction on k. 

When k = 0, it is clear that A(°), Q (0) , and i? (0) arc the same for both the simultaneous iteration method 
algorithm and the QR method algorithm, and these values satisfy ([1]) and J2j. 

Suppose now that the values of these matrices are the same for both algorithms for some iteration k — 1, 
and that they satisfy the two properties ([1]) and in this iteration. 

In the simultaneous iteration method, we have 

A k = AA k-i 

= A (Q**- 1 )^*- 1 )) 

= (q^r^)r^ 
= qWrW 

and thus property (JTJ) is satisfied on the fcth iteration; property ([2]) is satisfied directly by definition of the 
simultaneous iteration algorithm. 
In the QR method, we have 

A k = AA 1 "- 1 

= A (qWqWqW . . . Q**" 1 )^*- 1 ) . . . 

= (q^R^) (q^qWqW . . . Q(*-i)ij(*-i) . . . i?(D) 

= (Q( 2 >i?^) Q( 2 >Q( 3 > . . . Q(*-Dj2(*-i) . . . i?« 

= Q«Q( 2 > (q^ 3 ') Q< 3 > . . . qV>-»rV-i) . . . i?(D 

= Q^Q^V . . . Q {k -^ (Q«i?( fe )) R( k -V . . . 

= q^rW 
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which proves that property ([1} holds for the fcth iteration of the QR method. 
We also have that 

A W = Rik)Q{k) 

= ((Q (fc) ) T Q (fc) )i? (fe) Q (fc) 

= (qW) T 4qW 

which proves that property @ holds for the fcth iteration of the QR method as well. 

By hypothesis, the values of A^ fc_1 \ Q^ -1 **, and R^ 1 ' were the same for both algorithms. Since both 
algorithms satisfy (J]) on the fcth iteration, we have A k — Q^R} k > for both algorithms, and since the QR 
decomposition is unique, it follows that and are also the same for both algorithms. Both algorithms 

also satisfy ([2]) on the fcth iteration, which means that the matrix = {Q^ k ') T AQ^ k > is also the same for 
both algorithms on the fcth iteration. 

Hence, both algorithms produce the same values for the matrices A^ k \ Q , and B> k > which satisfy the 
relationships (JlJ and ((2|). By induction, this holds for all fc, proving the theorem. □ 

Now, since we saw that the columns of Q^ k ' converged to a basis of eigenvectors of A in the simultaneous 
iteration method, the result of the above theorem tells us that the same holds for the columns of Q^ k ' as 
computed using the QR algorithm. In particular, the columns q i of Q^ k ' each converge to a unit eigenvector 
qi of A. 

Property of the theorem tells us that 

where q^ and q^ are columns i and j, respectively, of Q^ k \ But, we saw that q^ — > qi and — > qj as 
fc — y oo , where qi and qj are unit eigenvectors of A and are orthogonal if i =fi j . 
In the case that i ^ j, we have that 

4 fc) -> qjAq, = X.qjq, = 

since the eigenvectors are orthogonal; thus A^ approaches a diagonal matrix. 
In the case that i = j, we have that 

A^ qf A Qi = X iq T qi = A, 

since the eigenvectors are unit vectors; thus the diagonal elements of A^> approach the eigenvalues of A. 

We finally see that the QR iteration method provides an orthogonal matrix Q^ k ' whose columns approach 
eigenvectors of A, and a diagonal matrix A^> whose diagonal elements approach the eigenvalues of A. Thus 
the QR method presents a simple, elegant algorithm for finding the eigenvectors and eigenvalues of a real, 
symmetric, full rank matrix. 

The QR method is, in essence, the same as the simultaneous iteration method, and therefore suffers the 
same restriction on the matrices it can be applied to. Unfortunately, these two methods may fail when 
there are non-real eigenvalues, or when the eigenvectors do not form an orthogonal basis of W l . However, 
symmetric matrices with real entries do occur frequently in many physical applications, so there is some use 
for these algorithms. 
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4. Conclusions and Improvements 

The purpose of this paper was to provide an overview of some of the numeric techniques used to compute 
eigenvalues and eigenvectors of matrices. These methods are all based on simple ideas that were steadily 
generalized and adapted until they became powerful iterative algorithms. 

However, even the QR iteration method as presented in this paper is generally not suitable for use in 
practice, even in the situations when it can be applied. There are many ways to refine the algorithms 
we have seen in order to speed up the implementations. For example, before applying the QR algorithm, 
it is common to reduce the matrix A into a simpler form, such as a tridiagonal matrix: performing QR 
decompositions on the new matrix then becomes much faster. It is also possible to improve the QR iteration 
method by incorporating shifts and Rayleigh quotients, just as these concepts helped improve the original 
power iteration method. The variations that are implemented in practice are vastly more sophisticated than 
the simple algorithms presented here. 

As with all algorithms, the question of improvement is always an open problem. Modern information 
processing deals with increasingly large matrices, and efficient methods for computing eigenvalues and eigen- 
vectors are extremely necessary. The techniques of the future may become much more complicated in order 
to keep up with this growing demand, but as we have seen, there is a lot that can be done using simple ideas 
and well-chosen generalizations. 
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